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Abstract. We report on the implementation of a coherent dipole shower algorithm along with an automated 
implementation for dipole subtraction and for performing powheg- and MCQNLO-type matching to next- 
to-leading order (NLO) calculations. Both programs are implemented as add-on modules to the event 
generator HerwigH — (-• A preliminary tune of parameters to data acquired at LEP, HERA and Drell-Yan 
pair production at the Tevatron has been performed, and we find an overall very good description which 
is slightly improved by the NLO matching. 

PACS. 12.38.Bx Perturbative QCD calculations - 12.38.Cy Summation of QCD perturbation theory 
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1 Introduction 

Many physics analyses at the Large Hadron Collider (LHC) 
are nowadays based on Monte Carlo simulations [IH5], 
e.g. for acceptance determination or even for background 
subtraction. With the high precision aimed for in many 
analyses it is mandatory to provide many of the simula- 
tions with the highest possible theoretical accuracy. For 
most processes this is now next-to-leading order (NLO) in 
the perturbative expansion of Quantum Chromodynamics 
(QCD). During the last decade, enormous progress was 
made in the development of techniques to match NLO 
calculations on the one hand and to merge multiple jet 
tree level matrix elements on the other hand with parton 
shower algorithms. 

First attempts to improve parton shower emission pat- 
terns with the information from the full matrix element 
for the hardest gluon emission were made with so-called 
matrix element corrections [BJE], that have long been im- 
plemented in the standard event generators. The next big 
improvement was made when matrix elements for multi- 
ple hard emissions were merged with parton shower al- 
gorithms, first for e + e~ annihilation processes [8j[9] and 
then also for hadronic collisions [ID] , An alternative ap- 
proach was proposed in |11] , where different implementa- 
tions have been systematically compared as well. The ex- 
perience that was made with these algorithms over the last 
years [H] has lead to further improvements [I31III] such 
that now the systematic uncertainties due to e.g. matching 
scale dependence have been significantly reduced. 

Matching to NLO matrix elements has been initiated 
first with a phase space slicing method [ToT - HT] . A more 
systematic matching has then been introduced by Frix- 
ione and Webber in the MC@NLO approach [TH]. This 
approach has then been generalised to include massive 



partons [19 . Many processes have been included in the 
meantime 20 22 . As the algorithm depends on subtrac- 
tion terms for a specific parton shower implementation, 
the first versions of MC@NLO have been tailored to work 
with Herwig only. Now, it also works with Herwig++, 
i.e. as the subtraction scheme has been generalised to- 
wards the HerwigH — I- parton shower implementation, all 
processes available in the MC@NLO package can also be 
showered with HerwigH — h to achieve formal accuracy at 

NLO Eg. 

As the matching of NLO matrix elements and parton 
shower algorithms takes place perturbatively to the speci- 
fied order, i.e. the next-to-leading order, there is formally 
an ambiguity left that can be used to devise alternative 
matching schemes. One such scheme has been proposed 
by Nason [51] and now goes under the name powheg. 
The guiding principle of this algorithm is to allow for a 
matching algorithm that does not introduce events with 
negative weight, as the MC@NLO prescription does. This 
approach has also been very successfully established dur- 
ing the last years and implemented as a separate program 
package [25] . Many processes are available in this program 
package [26-30 . However, the method itself is also used 
by other groups to match NLO calculations with parton 
showers within a given shower package. Many processes 
are available with Herwig+H- [3"TH3"5] or sherpa [55] , 
The internal implementations benefit from the inclusion 
of truncated showers (see below). 

On the parton shower side, a number of new parton 
shower algorithms have been developed during the last 
years, partly together with the rewrite of old generators 
37,38 . Many new developments have addressed the idea 
of implementing a shower that is directly related to the 
subtraction terms commonly used in NLO calculations. 
This led to the implementation of parton showers with 
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splitting kernels based on the Catani-Seymour subtrac- 
tion scheme g3HD] for NLO calculations gTJ|42], which 
was proposed in [33J . Similar ideas are followed with other 
subtraction schemes as e.g. in the VINCIA shower [33] where 
QCD antenna subtraction terms are facilitated. 

With more and more NLO calculations being matched 
one-by-one the question arises whether this step can be 
automated. In fact, the POWHEG method is already a first 
step into this direction, as the method as such is indepen- 
dent of the showering algorithm. In particular, no specific 
subtraction terms or the like are needed in order to match 
a given NLO calculation to any shower. There are sub- 
tleties on the shower side, though. The powheg method 
guarantees to give the hardest emission within the par- 
ton evolution and ensures that this is generated according 
to the phase space weighting of the NLO matrix element. 
However, if the shower does not evolve in the same hard- 
ness measure as the powheg algorithm, one has to intro- 
duce so-called truncated showers. This has been discussed 
already in early powheg implementations [35] and is now 
part of Herwig+-|- [H] and SHERPA [T3] . 

Many NLO calculations are available as ready-to-use 
computer codes that often come as packages that include a 
number of processes at NLO already. Most of these codes 
use the Catani-Seymour subtraction method to regularise 
infrared divergences. More recently, also the complete au- 
tomation of NLO calculations has been discussed with first 
tools readily available [351137] , based on the approach [48] . 
Some more calculations are already based on a fully au- 
tomated tool chain [4*9"h53"] . Part of this progress relies 
on the automatic generation of Catani-Seymour subtrac- 
tion terms [53H5S] or FKS subtraction terms [57]. The 
latest developments unify the matching of multiple tree- 
level emissions and the matching of NLO corrections to 
the Born level [581155]. 

In this paper we introduce an implementation of a 
parton shower based on the Catani-Seymour subtraction 
terms, similar to the showers introduced in [3UE2]- The 
goal of the implementation is to provide a framework for 
an automatic matching of NLO computations to a parton 
shower. The use of the subtraction terms is highly bene- 
ficial as the MC@NLO like matching, that is based on a 
subtraction of the parton shower contribution to the NLO 
observable becomes trivial. Together with a framework to 
handle powheg like matching we will have the possibility 
to check systematics within a single implementation. By 
using a shower based framework we may directly make use 
of truncated showers in order to minimise systematic un- 
certainties inherent to the matching formalism. As a first 
step in this programme we present the shower implemen- 
tation, which is embedded as a module in the HerwigH — h 
event generator. In addition we present NLO matchings to 
the basic QCD processes. 

The paper is organised as follows. In Sec. [5] we intro- 
duce the dipole shower in detail. Sec. [3J introduces the 
implementation of an automatic matching with this par- 
ton shower, that we call Matchbox. In Sees. [4] [5] and |6j 
we present comparisons to data from LEP, HERA and the 
Tevatron, respectively. 



2 Dipole Showers 

The dipole shower algorithm outlined in [60] has been im- 
plemented as an add-on module to HerwigH — h, PQ. In 
this section we briefly review its properties and give a full 
description of the implementation. 

The authors have shown that parton showers based 
on Catani-Seymour subtraction kernels [39] correctly re- 
produce the Sudakov anomalous dimensions and properly 
include effects of soft gluon coherence, upon using an or- 
dering of emissions in transverse momenta as defined by 
the emitting dipoles. The simple inversion of the kinematic 
parametrisation used in the context of NLO subtraction, 
however, does not resemble a physical picture for initial 
state radiation. An alternative has been suggested and 
implemented in the simulation presented here. 

2.1 Starting the Shower 

The dipole shower starts evolving off a hard sub process, 
which is assigned colour flow information in the large- N c 
limit. This colour flow information is used to first sort all 
coloured partons attached to the hard sub process into 
colour singlets. Practically, this is done by making use of 
the fact that a colour singlet is 'simply connected' in the 
sense of its colour flow topology: Any parton i in a colour 
singlet can be reached from a parton j in the same singlet 
by just following colour lines and changing from a colour 
to an anti-colour line at an external gluon. Each colour 
singlet is now an independently evolving entity, and can 
only split into two colour singlets in the presence of a 
g —} qq splitting. In the next step, the partons in each 
singlet are sorted such that colour connected partons are 
located at neighbouring positions, when representing the 
singlet group of partons as a sequence. Note that these 
sequences may be open or closed: We will call a sequence 
open, or non-circular, if there exists a circular permutation 
of the elements in it such that the partons at the first and 
last position are not colour connected. Conversely, if there 
does not exist such a permutation, the sequence is called 
circular or closed. The possible sequences are depicted in 
Fig. [1] Once this sorting has been accomplished, we will 
refer to these singlet sequences as dipole chains: each pair 
of subsequent partons in a singlet sequence forms a dipole, 
which may radiate. For each parton in each dipole, a hard 
scale is then determined as defined in |60j . 

2.2 Evolution of the Parton Ensemble 

The main shower algorithm acts on a set of dipole chains, 
and proceeds as long as this set is not empty. Dipole 
chains are removed from the list, if they stopped evolv- 
ing, i. e. if there was no splitting selected with a p\ above 
the shower's infrared cutoff ^j R . The first entry in the 
set of dipole chains is taken to be the current chain. For 
each dipole in the current chain (with both possible 
emitter-spectator assignments, i.e. also considering (J,i) 
along with any possible splitting (i,j) — > (i',k,j) 
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is considered to compete with all other possible splittings 
of the chain. For any such splitting, given a hard scale 
p\ associated to the emitter under consideration, a scale 
q\ is selected with probability given by the Sudakov form 
factor 

exp ( - / dq 2 I dzP( i!i) _ ) .( i / ifcii )(3 2 ,z) , (1) 



where Pi 



-(i 2 ) 



(q 2 , z) is the appropriate splitting prob- 



ability as defined in [5D] , using the respective dipole split- 
ting function V5',fe;i< 

The splitting with the largest selected value of q\ is 
then chosen to be the one to happen, except the largest 
q\ turned out to be below the infrared cutoff. In this case 
the current chain is removed from the set of dipole chains, 
inserted into the event record and the algorithm proceeds 
with the next chain. The momentum fraction z is chosen 
to be distributed according to dP^ ,j)-*(i',k .j) (<0_) 2 0- Since 
for now we use azimuthally averaged splitting kernels, the 
azimuthal orientation of the transverse momentum is cho- 
sen to be distributed flat. The momenta of the splitting 
products and the spectator after emission are then calcu- 
lated as specified in [60] . 

As the evolution factors into dipole chains as indepen- 
dently evolving objects, all possible emitters in the chain 
- after having inserted the generated splitting - now get 
the selected q\ assigned as their hard scale, or stay at 
the kinematically allowed scale p\ i ■ if q\ > p\ i ■ . If a 
g — > qq splitting has been selected for a circular chain, this 
chain becomes non-circular. If it has been selected for an 
already non-circular chain, this chain breaks up into two 
independent chains exactly between the gg-pair, owing to 
the colour structure of this splitting. This situation, along 
with non-exceptional splittings is depicted in Fig. [1] 



2.3 Finishing the Shower 

After the shower evolution has terminated, the incoming 
partons with momenta p a ^ in general have non- vanishing 
transverse momenta with respect to the beam directions. 
This necessitates a realignment of the complete event en- 
countered at this stage. Following the arguments of [60 , 
the momenta of the evolved incoming partons p a .b are 
taken to define the frame of the collision at hand, i.e. 
hadron momenta P a ,b- We then seek a Lorentz transfor- 
mation to take Pa t, to the externally fixed hadron mo- 
menta P a .b, which is in turn used to realign the complete 
event. 

To construct the momenta of the incoming hadrons 
Pa,bi we require the three-momenta of P a ,b being collinear 
with the respective partonic three-momenta and define 
momentum fractions 



•Ea.b 



2-Pfe.o • Pa.b 



o 



o. 



Gluon emission off a circular chain. 
The chain stays circular. 



o. 



Gluon emission off a non-circular chain. 
The chain stays non-circular. 



o 



Q 



a 



g — > qq splitting in a circular chain. 
The chain becomes non-circular. 



Q 



^1 

O 



o: 



(2) 



g — > qq splitting in a non-circular chain, 
triggering breakup of the chain. 

Fig. 1. Examples of parton emission from dipole chains. In 
these examples always the upper dipole has been considered 
for emissions. Note that any dipole may split in two different 
ways, splitting either of its legs. These competing possibilities 
are not shown in the transition diagrams. 



The momentum fractions are further constrained by re- 
quiring that 

(Pa + Pb) 2 = S (3) 

where S is the centre-of-mass energy squared of the colli- 
sion, such that the desired Lorentz transformation exists. 

The second constraint is in principle to be chosen in 
such a way as to preserve the most relevant kinematic 
quantity of the hard process which initiated the showering. 
By default, we choose this to be the rapidity of a system 
X , which is cither the system of non-coloured particles at 
the hard sub-process, or the complete final state in case 
of a pure QCD hard scattering. 



2.4 Cluster Hadronization 

The cluster hadronization model, originally proposed in 
[6"T] . is the hadronization model used by the Herwig+ + 
event generator. The model in its initial stage just af- 
ter parton showering, performs a splitting of gluons into 
quark-antiquark pairs such that in the large-A^ limit a 
set of colour singlet clusters emerge from the event under 
consideration. 

These clusters are then subsequently converted into 
hadrons, by either splitting them into clusters of lower 
invariant mass or performing directly the decay to meson 
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pairs, in case another qq pair is 'popped' from the vacuum 
inside the cluster, or baryon pairs, where the creation of 
a diquark-antidiquark pair is assumed. Further details of 
the model will not be discussed here. 

The main assumption of the model is however, that 
both quarks are located on their constituent mass shell, 
and gluons are as well assigned a non- vanishing constituent 
mass, entering as a parameter of the model. In the stan- 
dard HerwigH — h parton shower, acting as a 1 — > 2 cas- 
cade, only scales and momentum fractions of the splittings 
are determined during the evolution, the full kinematic 
information being constructed after the end of the pertur- 
bative evolution. This setup thus straightforwardly allows 
to include the constituent masses in this particular step. 
Since the dipole shower preserves momentum conservation 
locally to each splitting, ending up with a set of massless 
partons, such a treatment is not possible. 

The way to perform the 'reshuffling' of the massless 
parton momenta to their constituent mass shells is chosen 
to be the following algorithm: Let Q c be the total momen- 
tum of all final state partons and perform a boost A c to the 
centre-of-mass system of Q c , A C Q C = (Q c , 0). The boosted 
parton momenta pi are now put on the constituent mass 
shell, including a global rescaling of their three-momenta, 



Pi 



(IPil.Pi) ^ri = (^ 2 |Pil 2 + "£i,£p. 



(4) 



Momentum conservation then implies the following rela- 
tion be satisfied, 



E 



(5) 



which may be solved numerically to yield a value for £. 
Finally the inverse boost A' 1 is applied to the new parton 
momenta p\. 



3 The Matchbox Implementation 

Closely related to the dipole shower implementation, though 
technically independent of it, is the development of the 
Matchbox module. Matchbox is based on an extended 
version of ThePEG, the extensions providing functional- 
ity to perform hard process generation at the level of NLO 
QCD accuracy and easing the setup of run time interfaces 
to external codes for hard process generation. We have im- 
plemented an automated generation of subtraction terms 
based on the dipole subtraction formalism [39], based on 
the information available from ThePEG matrix element 
implementations, which will be discussed in further de- 
tail in section 13.21 A full NLO calculation to be run in 
the Matchbox framework only requires the implementa- 
tion of tree-level and one-loop amplitudes, the presence of 
colour (and spin) correlated amplitudes for the Born pro- 
cess and the presence of a phase space generator appro- 
priate to the process under consideration. Fig. [2] sketches 
the involved software modules and their interaction with 
an external implementation of a NLO calculation. 



NLO Code 
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Hadronization and 
Underlying Event 



Dipole Shower 



NLO Setup 



Matchbox Dipoles 



Matching Setup 



ExSample 



Fig. 2. A sketch of the interaction of the Matchbox and dipole 
shower modules as integrated in HerwigH — h To perform a 
matched NLO calculation an external code only has to provide 
tree-level and one-loop amplitudes along with colour- and spin- 
correlated amplitudes of the Born process and an appropriate 
phase space generator. 



Besides being capable of performing a Monte Carlo in- 
tegration of 'plain' NLO corrections, the main purpose of 
Matchbox is to turn a NLO calculation into a matched 
calculation to be consistently combined with a parton 
shower. Here, functionality is especially provided to cal- 
culate the inclusive NLO cross section differential in the 
Born degrees of freedom, which, along with a matrix ele- 
ment correction to the shower, is the main ingredient to 
the POWHEG method of combining parton showers and 
NLO QCD corrections. 

Matchbox is automatically generating matrix element 
corrections from the NLO real emission contribution. It 
further allows the possibility to overcome problems in the 
powheg matching owing to radiation zeroes in the Born 
matrix element. The matrix element correction splitting 
kernel, which is essentially defined by the ratio of real 
emission and Born matrix elements squared is turned into 
the corresponding distribution including the Sudakov form 
factor by using the ExSample library, [52]. ExSample 
allows the efficient sampling of distributions of this type, 
without having to provide any analytic knowledge on the 
splitting kernel or trying to estimate enhancement fac- 
tors to simpler functions such as dipole splitting kernels. 
ExSample is also used to sample emissions in the dipole 
shower implementation. 

3.1 Notation 

We consider NLO calculations carried out using the dipole 
subtraction method, [32] • Instead of using the notation es- 
tablished there, we unify the indices of all possible dipoles 
to ease readability, as expressions become quite compli- 
cated especially when considering the powheg type match- 
ing. For the subtraction dipoles we choose the notation 



V. 



ij,k 



V a - 

> v 



ai,b 



— > 



(G) 



where the arguments are unified and we make explicit the 
dependence on either real emission or 'tilde' kinematics, 
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e.g. 

T^ij,k(q a ,qb] qi, ■••) q n +i) 



T> a {p"(q n +i)\q n +i) ■ (7) 



In this notation, p n now refers to the whole phase space 
point, 

p a ,p b ;p 1 ,...,p n -> p n = (p a ,p b ;p 1 ,...,p n ), (8) 

where we have added hat symbols to the momenta to dis- 
tinguish a single momentum from a complete phase space 
point. The 'tilde' mapping and its inverse are denoted by 

Pij{q l ,qj,qk) ,Pk(qi,qj,qk) -*v>%{q n +i) (9) 

q l .j,k(p7j,Pk;p 2 L ,z,(t)) -> q" +1 (p n ;p 2 ± ,z,(j)) . 

Differential cross sections are considered in collinear fac- 
torisation, 

d<T X (Pn\Q,Xa,X b ,H F ) = 

fp^a{x a , HF)fp^b{xb, HF)d<7 X (Pn\Q)dx a dx b (10) 

where the partonic cross section is in general of the form 

da x (Pn\Q) = F{p a ,p b )X(p n )d<fi(p n \Q) . (11) 

Here F(p a ,p b ) is the appropriate flux factor and X(p n ) 
generically denotes any contribution to the cross section 
which can be cast in the above form, i. e. tree- level ampli- 
tudes squared, one-loop tree-level interferences, subtrac- 
tion terms, or the 'deconvoluted' finite collinear terms to 
be discussed below. The phase space measure d<fi(p n \Q) is 
given by 



d4>{Pn\Q) 



(ibr)'* (J>~Pa -P„ ~ Qj g ^L^o (12) 

In latter sections, it will turn out to be useful to rewrite 
this as 

da x (Pn\Q,x a ,x b ) = X(p n )dF(x a ,p a ,x bl p b )d(j){p n \Q) 

= X(p n )d(j) F {p n \Q,x a ,x b ) , (13) 

where we dropped making explicit the factorisation scale 
dependence from now on. 

The finite collinear terms originating from counter terms 
to renormalise parton distribution functions and integrated 
subtraction terms are reported in |39j . These are given as 
convolutions of Born-type cross sections of colour corre- 
lated amplitudes with certain 'insertion operators', e.g. for 
the incoming parton a 



events should be generated according to the rescaled in- 
coming momentum zp a . A numerical implementation is at 
first sight not obvious. Considering however the integra- 
tion over the momentum fraction x a , these contributions 
can be rewritten in terms of a Born-type cross section 
multiplied by modified PDFs along the lines of 



dx [ dzf(x)B(xz)P(z) 



dz C(p a n {z))d(j>{p n \Q a (z))dF(x a ,zp ai x b ,p b ) , (14) 



where the superscript a along with an argument z indi- 
cates, that parton a's momentum is rescaled by z. The in- 
sertion operators themselves include -(--distributions, and 







\ xB {x) r^f^p( Z ) (is) 



and the +- distributions can be expressed in a way to allow 
for numerical implementation. All possible contributions 
for light quarks are implemented in Matchbox. 

Any NLO cross section within the dipole subtraction 
thus takes the form 

vnlo = J \M B (p n )\ 2 u{p n )d<j) F {p n \Q,x a ,x b ) (16) 

+ J [2R C {M* B {Pn)Mv(Pn)) + 

{M B {Pn)\I\M{p n ))] e=0 u{p n )d<t> F {p 

+ / (M B (Pn)\(P + K)\M( Pn ))u( Pn )dMP 

+ J (\M R {q n+ i)\ 2 u{q n+1 ) 

- ^2v a (p%(q n+ i)\q n+1 )u(p%(q n+1 )) J 

a / 

x d(/) F (q n+ i\Q,x a ,x b ) 

where the insertion operators I are given in [39J and have 
been implemented for light quarks in full generality as 
well. P, K and d<p F denote the deconvoluted versions of 
the finite collinear terms originating from the insertion op- 
erators P,K given in [35]. Here, the test functions u{p n ) 
refer to the class of events to be generated by a Monte 
Carlo realisation of the above integrals, and M b ,r de- 
note the Born and real emission amplitudes, respectively. 
Since only the structure of the real emission and subtrac- 
tion terms turns out to be relevant for matching purposes, 
we from now on collectively denote Born, virtual and in- 
sertion operator contributions by 

\M B v(pn)\ 2 u(p n )dci) F (p n \Q,x a ,x b ) . 

Since all the integrals will be dealt with by means of 
Monte Carlo methods, differentials are expressed in terms 
of a Jacobian expressing the physical variables in terms 
of random numbers and a volume element on the unit 
hypercube of these random numbers, e.g. 



d<t>(p n \Q) 



dp„ 



Or 



d k r 



(17) 



We identify ratios of differentials to actually mean the ra- 
tios of the corresponding functions multiplied by the Jaco- 
bian in use to express them in terms of random numbers, 
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e.g. for two cross sections we define 



d(T X (q m \Q) 



da Y (p n \Q) 



X(q m ) 


dq m 


dr q 


Y( Pn ) 


dp n 
dr p 



(18) 



3.2 Automated Dipole Subtraction 

Any matrix element implemented in ThePEG is expected 
to provide information on the diagrams contributing to 
it. It is this information, which is used to generate sub- 
traction dipoles by a simple algorithm of checking, for any 
contributing diagram, if any two external coloured legs are 
attached to the same vertex. By removing this vertex from 
the diagram information, the diagram of the correspond- 
ing 'underlying Born process' is obtained. Conversely, the 
same pairing of diagrams provides a way to identify which 
real emission processes are to be considered given any 
Born process. This information is used when setting up 
the inclusive NLO cross section calculation and generating 
matrix element corrections for the parton shower. From a 
given matrix element object implementing a real emission 
contribution, Matchbox checks a set of Born matrix el- 
ement objects provided along with the real emission ones 
for the underlying Born processes obtained and adds all 
matching pairs to the calculation if there exists a sub- 
traction dipole object which claims responsibility for the 
given pairing. Similarly, all insertion operator implemen- 
tations present are checked if they claim responsibility for 
a given Born process, thus completing the setup of a NLO 
calculation. The complete calculation is then injected as 
a ThePEG SubProcessHandler object into the stage of 
event generation. 

For running unmatched calculations, a group of events 
consisting of real emission and 'tilde' phase space points 
is provided along with the relative weights of the individ- 
ual contributions present in the group. The sum of these 
weights, i.e. real emission minus subtraction term contri- 
butions is driving the cross section integration and event 
unweighting. 



3.3 Subtractive NLO Matching 

Owing to the fact that the dipole shower implementation 
uses splitting kernels which precisely equal the dipole sub- 
traction terms, following the steps leading to MC@NLO 
here results in a very simple matching^ This subtractive 
matching is basically identical to the NLO calculation it- 
self, except that instead of event groups now a single real 
emission phase space point is generated from the sub- 
tracted real emission contribution. In an algorithmic man- 
ner, the matching may thus be expressed very simply: 



1 Though the kinematic parametrisation differs from the one 
used in the subtraction context, it can be related to the usual 
'tilde' parametrisation by a boost in case a single emission is 
considered. 



— Generate Born-type events p n with density 

\MBv{Pn)\ 2 A(j)F{Pn\Q,X a ,X b ) , (19) 

— generate real-emission type events q n +i with density 

(\M R (q n+1 ) | 2 - £ V « (<Zn+l) k+l)J 

x d(f> F {q n+1 \Q,x a ,x b ) , (20) 

— and feed either into the dipole shower. 

A subtlety, however, arises here. Since we are interested 
in describing the hardest emission according to the exact 
real emission matrix element, the parton shower should 
not generate harder emissions than the one fixed from the 
NLO calculation. Practically, this is implemented by cal- 
culating the p°l_ as defined by the inverse 'tilde' mapping 
from each dipole configuration a, since the kinematics of 
the emission appears differently depending on the emitting 
dipole considered, is communicated as a veto scale to 
the dipole shower, which is not allowed to generate emis- 
sions with p± > pj_ off the emitter, emission and spec- 
tator partons used to evaluate T> a . Another approach, in 
which the dipole shower is generally not allowed to emit 
at scales p±_ larger than final state transverse momenta 
can equivalently be used and may become the default in 
a future version. This treatment is then very similar to 
the Herwig shower in use with the traditional MC@NLO 
implementation. 

3.4 NLO Matching with Matrix Element Corrections 

The splitting kernels to be used for a matrix element cor- 
rection are given by the ratio of real emission and Born 
matrix elements squared, weighted by (in principle) ar- 
bitrary weight functions for each kinematic mapping of 
a subtraction term, i.e. for each subtraction term. It is 
most simple to choose the subtraction terms themselves 
to define these weight functions. This has the advantage 
that all divergences but the divergence associated to the 
subtraction term D a are divided out from the real emis- 
sion matrix element, and dynamical features of the Born 
matrix element, like peaks owing to unstable particles, are 
flattened out in the splitting kernel considered. 

Within this procedure, one faces three major problems: 

— Some of the subtraction dipoles, in particular the ones 
with initial state emitter and final state spectator or 
vice versa, are not positive-definite. This makes a Monte 
Carlo treatment of the corresponding Sudakov-type 
distribution hard to implement. Since the regions, where 
these dipole kernels become negative correspond to 
hard, large angle parton emission, it is clear that this 
problem can be cured by changing the irrelevant fi- 
nite terms of the subtraction dipoles, provided they 
are consistently taken into account in the integrated 
ones. Within the Matchbox implementation this has 
so far been carried out for the qq initial-final dipoles, 
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which have been modified to reproduce the matrix ele- 
ment squared for gluon emission off the corresponding 
vector current and are thus positive by definition. 

— The Born matrix element squared may contain zeroes. 
In this case, its inverse is obviously ill-defined. 

— The implementation of the parton densities at hand, 
which enter as a ratio in the splitting kernels as well, 
may not be stable in particular for large x in the sense 
that the interpolation used oscillates around zero rather 
than tending to zero smoothly. This poses a problem 
similar to the zeroes in the Born matrix element, how- 
ever now without any physical interpretation. 

The latter two problems can be solved by introducing 
an auxiliary cross section da SC reen(pn\Q', P±) which enters 
into the definition of the splitting kernels 



and 



dPa(p 2 A_,z,4>\Pn) = d 3 r 



E^(K(C+iM+i) 

da R (q% +1 \Q 



d<7 B {Pn\Q ) +da i 



screen, a 



(Pn\Q;p±) 



(21) 



where we have already written the splitting kernel differ- 
ential in the random numbers determining pj_, z and <fi, 
and the dependence of q^+i = 9n+i(PniPj_> z i 4>) on the 
splitting variables is understood implicitly. In order not 
to change the divergence structure implying the resumma- 
tion of large logarithms, the screening cross section needs 
to vanish as p\ — > 0. Since Born zeroes cannot occur for 
p\ — >• (the QCD singularities factor in this limit with 
respect to the Born process) Eq. (l2~Tj) is free of these prob- 
lems. If, in addition, the screening cross section does not 
depend on the parton distributions, the technical issues 
with PDFs becoming zero are cured as well. 

The screening cross section has however to be taken 
into account for the fixed order calculation in order to 
reproduce the correct NLO cross section and will thereby 
spoil the original simplicity of using the NLO if-factor dif- 
ferential in the Born variables to generate events to enter 
the matrix element corrected shower. Including the screen- 
ing cross section the fixed order cross section can then 
be calculated to be constructed of densities for Born-type 
and real emission type events. The densities for Born-type 
events closely resemble the X-factor modification, 

d(Tinclusive {jpn | Q : X a , Xfr) — 

i / l/^i \ , f j3 dtXR inclusive {Pn \ Q> %ai x b) 

da B v{Pn\Q,x a ,x b ) + / d r 



d 3 r 



(22) 



where 

d<7 R, inclusive (.Pri \ Q j %a j %b ) 



d fe 



d k r B d 3 r 



da B (p n \Q) Pq(P«|gg + i) 



(23) 



RiPnK+x) 



d(f>(p n \Q) 
da R (q% +1 \Q 



da B (Pn\Q,X a ,X b ) + d(T scre en,a(Pn|Q;Pi) 



(24) 



To generate events according to these densities, a k + 3- 
dimensional random number point is chosen, where the 
three additional degrees of freedom are discarded. Owing 
to the fact that the integration volume in terms of random 
numbers is the unit hypercube, this procedure produces 
the integration over the degrees of freedom of the parton 
emitted in the real emission on average. 

Events of real emission type are to be generated with 
density 

da R {q n+1 \Q,x a ,x b ) x 

R(p n \q n+1 ) , (25) 

2^0 V l3{Pn\<ln+l) 



R(Pn\ln+l) 



da*. 



da B (p%\Q (Pn\Q;p±) 



(26) 



which is just a reweighting of the real emission contribu- 
tion. Events of both classes can then be showered by a 
parton shower using a matrix element correction as de- 
fined at the beginning of this section, and a communica- 
tion of veto scales applies to the real emission contribution 
along the same lines as for the subtractive matching. Note 
that the individual contributions are positive, as long as 
the screening cross section is bounded from above by a 
reasonable value. 

Since this type of matching is independent of the par- 
ton shower to act downstream, the actual implementation 
does not make any reference to the dipole parton shower, 
and real emission contributions according to the matrix 
element correction are generated outside any shower mod- 
ule, presenting a real emission sub process supplemented 
with proper veto scales, or a Born-type sub process to the 
shower, if radiation has been generated according to the 
matrix clement correction or not, respectively. 

Note that, when putting the screening cross section to 
zero, the original simplicity of the POWHEG-type match- 
ing is recovered. The matrix element corrections, inclusive 
and real-emission type contributions are all setup and cal- 
culated in an automated way within the Matchbox im- 
plementation. The screening cross section is by default 
chosen from the corresponding phase space and the di- 
mensionality required by the phase space, i.e. 



der Sl 



,, a (Pn(<ln+l)\Q;pl) = 



(Pa 



d(f>(q n +i\Q) 



t (q n +i) (s a (q n+ i)) n ™ 



(27) 

where is the transverse momentum associated to the 
mapping (q n +i ) , s a {q n +i) is the appropriate mass squared 
of the emitter-spectator pair in p%, and n out is the number 
of outgoing particles. Other choices may be possible. 
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4 Results at LEP 

The variety of data acquired by the LEP experiments al- 
low for a systematic fit of parameters of the parton shower 
and the hadronization model. In a preliminary fit, the pa- 
rameters assumed to mainly determine the description of 
event shape variables and jet rates as measured by the 
DELPHI experiment [53J and jet observables as reported 
by the OPAL collaboration [51] have been fitted using the 
Rivet [55] and Professor [55] systems. The parameters 
and ranges considered are given in Tab. [TJ along with a 
short description. Parameters which are known to mainly 
affect individual hadron multiplicities have not been var- 
ied, and fragmentation parameters for heavy quarks have 
been set equal to the values of those for light quarks. A 
simple modification of the running of a s in the infrared has 
been adopted by replacing its argument q 2 — s> q 2 + ^ 2 oft - 
This modification has originally been motivated to supply 
another model for intrinsic transverse momentum gener- 
ation by letting the initial state shower evolve down to 
very small scales along the lines of [57]. We see however 
no reason that it should not be considered for final state 
radiation as well. 

Separate fits have been performed for LO and NLO 
predictions. LO predictions have been obtained by run- 
ning just the parton shower, using a one-loop running a s . 
NLO prediction have been obtained by means of supple- 
menting the shower with the matrix element correction 
matching without using the Born screening cross section 
and a two- loop running a s . In total we find that the NLO 
simulation gives a marginally better fit than the LO one, 
though the description of data is completely comparable 
within experimental uncertainties. 

The fitted parameter values are displayed in Tab. [5] 
Most notably, the hadronization parameters for the LO 
and NLO fit do not significantly differ. For both predic- 
tions, a modification of the infrared running of a s seems 
not to be preferred. The infrared cutoff of the parton 
shower is determined more precisely by the NLO fit, which 
prefers a smaller cutoff. Also a s (M\) is determined more 
precisely by the NLO fit. Both a s values obtained are 
compatible with the world average [68 of 0.1184, where 
the NLO result is closer to this value. Note that this 
should be regarded a coincidence at the level of the ap- 
proximation considered and it is certainly not possible to 
uniquely relate the obtained value to one applying to the 
MS scheme. In Figs. [3J and H the LO and NLO simula- 
tion results are compared for selected observables. Fig. [5] 
shows the energy-energy-correlation, which has not been 
included in the fit. 



4.1 Comparison of Matching Strategies 

The Matchbox framework provides the facility to switch 
between the POWHEG-type matching with matrix element 
corrections including or excluding the auxiliary Born screen- 
ing cross section, and subtractive matching. For reasons of 
systematics it is instructive to compare these approaches. 
No separate fit for the variants not considered so far has 
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Fig. 3. Some event shape variables as predicted by the leading 
order and next-to-leading order simulations. Here, we addition- 
ally compare to the standard HerwigH — h shower (version 2.5.1 
with default settings), showing that the dipole shower gives a 
significantly improved description already at leading order. 



Simon Platzer, Stefan Gieseke: Dipole Showers and Automated NLO Matching in Herwig++ 



9 




Energy-energy correlation, EEC 



DELPHI data 
Dipole Shower 
Dipole Shower + NLO 




Fig. 5. Energy-energy correlation. Note that this observable 
has not been included in the fit. 



been performed and the NLO fit values as given in the 
previous section have been used. The different matching 
strategies give completely comparable results. If there are 
small visible differences, there is no clear tendency that 
either variant would give a better description than any 
of the others. Fig. [6] compares the matching strategies for 
the two jet rate. To this extent, the subtractive matching 
could be preferred amongst the POWHEG-type ones owing 
to its smaller computational complexity. This statement, 
however, not only includes that negative weighted events 
do not pose a major problem, but also has to be verified 
in a process dependent matter since there is no hint, if the 
behaviour observed here is a general feature - particularly 
at hadron colliders. 



a given dipole may be considered a gauge invariant 
subset in the soft and/or collinear limits for N c — ^ oo. 
This implies that the infrared cutoffs and soft scales en- 
tering the emission probabilities need not be the same for 
all dipoles. The emitter-spectator configurations forming 
gauge invariant quantities in this sense are the two emit- 
ter choices for final-final dipoles, initial-initial dipoles, and 
the combination of initial-final and final-initial configu- 
rations. Fitting DIS data therefore allows one to fix the 
infrared cutoff and soft scale for the latter, before finally 
constraining the same parameters for initial-initial dipoles 
at a hadron collider, which is considered in the next sec- 
tion. 

For the fit described here, the same technique as for 
LEP, and data accumulated by the HI experiment [69. 
have been used. For LO and NLO, the default Herwig++ 
PDFs, MSTW 2008 LO** [7DH7T] and MRST 2002 NLO 
[72"] . have been used. The same PDFs were considered for 
hadron collider data to be discussed in the next section. 
The NLO fit was obtained by running the matching with 
matrix element correction. 

The findings are similar as for the fit to LEP data. 
We find a reasonable prediction of transverse energy flows 
over the whole range of (x, Q 2 ) plane. The matched NLO 
prediction gives a comparable fit to the LO simulation, 
while preferring both a smaller infrared cutoff and screen- 
ing scale. The fitted parameters are given in Tab. [3] 

Fig. [7J shows the average transverse energy as a func- 
tion of Q 2 in the central detector region. This observable 
is clearly improved by the NLO matching at small mo- 
mentum transfers. A more detailed analysis of DIS data 
including inclusive jet and event shape data is currently 
underway. 



5 Results at HERA 

Owing to the approximation underlying the dipole par- 
ton shower, diagrams contributing to parton emission of 



6 Results at the Tevatron 

After having determined the simulation parameters for 
hadronization, final state radiation, and radiation off a 
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Fig. 7. Average transverse energy in the central region as mea- 
sured at HERA and compared to leading order and next-to- 
leading order predictions. 



Fig. 8. Differential cross section of the Drell-Yan-pair p± com- 
pared to LO and NLO predictions. Note that the cross sections 
have been normalised to the measured one. 



Pscudorapidity, of 3rd jet 



final-initial dipole by fitting LEP and HERA data, two pa- 
rameters remain to be determined: the infrared cutoff and 
soft scale for radiation off an initial-initial dipole. We here 
consider the p± spectrum of e + e~ Drell-Yan pair produc- 
tion as measured by the CDF collaboration [73] . Since the 
Drell-Yan process receives rather large QCD corrections 
from leading to next-to-leading order and a still consider- 
able correction at NNLO, both fits have been performed 
by normalising the simulation to the measured cross sec- 
tion. The matrix element matching including the Born 
screening cross section has been used here, as for the DIS 
data. 

The Professor algorithm here turned out not to be 
applicable, as the cubic interpolation was not capable of 
describing the complete dynamics of letting the shower 
evolve to rather small infrared cutoffs, owing to the pre- 
scription of introducing a soft scale in a s as already de- 
scribed before. We have therefore performed a prelimi- 
nary fit by generating 300 random points uniformly in 
parameter space, which here includes the infrared cutoff 
for initial-initial dipoles, the soft scale for initial-initial 
dipoles, as well as the widths of a Gaussian distribution 
for intrinsic transverse momentum, A±. The latter has 
been chosen to be potentially different for valence and sea 
partons. 

Out of these random points we have picked the one 
with lowest x 2 with respect to the data, again both for LO 
and NLO simulations. The resulting parameters are given 
in Tab. |H Note that the p± distribution for sea partons is 
narrower, corresponding to a broader spatial distribution 
as can be motivated on different grounds. 

We show the comparison of LO and NLO simulations 
in Fig. [5] showing similar systematics to the distributions 
discussed before. In order to determine the predictivity of 
the simulation already at this very coarse level of tuning, 
we additionally show the pseudo-rapidity distribution of 
a third jet in events with at least two hard jets, Fig. |H1 as 
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Fig. 9. The pseudo-rapidity distribution of a third jet in 
events with at least two jets. We here only show the leading 
order prediction in order to check the predictivity of the tune 
carried out so far. 



carried out at CDF [71] . Reasonable agreement with data 
is found. On top of the work presented in [50], this con- 
stitutes another crucial test of coherent parton evolution. 



7 Conclusions 

We have introduced a new dipole shower module for the 
event generator HerwigH — \- that allows for an automatic 
matching of NLO computations with a parton shower. 
A tune of the hadronization module to the most impor- 
tant data sets show that we can achieve very good results 
from this simulation already without the inclusion of NLO 
terms. Including NLO corrections at this relatively simple 
level only marginally improves the results. This effect is 
expected as it is known that the Catani-Seymour showers 



Simon Platzer, Stefan Gieseke: Dipole Showers and Automated NLO Matching in Herwig++ 



11 



tend to mimic the behaviour of NLO matrix elements very 
well also in phase space regions well outside the collinear 
limits. However, the matching poses no technical prob- 
lem and can be seen as a proof-of-concept for the idea 
to provide a framework for automatic matching. At this 
time with relatively simple matrix elements at NLO that 
are provided by internal code. Future work will concen- 
trate in the inclusion of external code via a well defined 
interface, following the ideas in [75] . 
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A Code Validation 

A.l Shower Splitting Kernels 

The sampling of shower splitting kernels has been explic- 
itly verified in situ, meaning using the full implementation 
as present in the simulation code, against an independent 
implementation using a numerical integration to obtain 
the Sudakov-type distributions. Fig. [TU] shows an example 
for a final-final splitting kernel, proving correctness of this 
part of the code. 

A.2 NLO QCD Corrections 
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Fig. 10. Example comparison of sampled final-final splitting 
momentum fraction (blue lines) versus results from a numeri- 
cal integration (turquoise lines) at two different dipole masses, 
Sij = (lOOGeV) 2 (continuous lines) and Sij = (50GeV) 2 (bro- 
ken lines). 
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All leading order matrix elements implemented in the MATCHFig. 11. Envelopes of the ratio of the subtraction to the real 
BOX framework have been cross-checked against the Her- emission cross section versus the propagator denominator for 
WIG++ matrix elements. all singular configurations in Z + jet production. 

The functionality of the automatically generated sub- 
traction terms has been verified. Fig. [TT] shows a typical 
examples of the ratio of subtraction to real emission cross 
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section, plotted against each of the invariants entering the 
propagator denominators. 

The 'plain' NLO cross section, and the inclusive one 
entering the matching with matrix element correction have 
been checked to agree, with and without the usage of the 
Born 'screening' cross section. The NLO cross section for 
e + e~ — > jets has been validated against the analytically 
known X-factor of 1 + a s /iv. The NLO cross section for 
DIS and Drell-Yan has been checked against the existing 
powheg implementation in Herwig+ + - For deep inelas- 
tic scattering, the subtraction terms have been modified in 
order to have positive definite dipole kernels, finite terms 
of the integrated subtraction terms have been changed 
accordingly. The functionality of the subtraction has been 
checked with both variants, and the NLO cross sections 
with and without modifications are found to agree. 



A. 3 NLO Matching with Matrix Element Corrections 

A non-trivial cross check of the matrix element correction 
code and ExSample as the underlying 'working horse', is 
to consider the spectra for a gluon emission off a qq dipole 
as generated by the shower, which is validated against 
a numerical integration of the expected distribution im- 
plemented in a completely independent code. By putting 
the real emission matrix element entering the matching to 
be equal to the sum of dipoles (the correctness of which 
has been checked by verifying that the cross section of 
the subtracted real emission matrix element is consistent 
with zero), the matrix element correction must produce 
the same spectrum as the shower code. We have checked 
that this is indeed the case. It should be stressed that the 
machinery underlying the setup of the matrix element cor- 
rection is much more complex than the shower implemen- 
tation, and, that the splitting kernel entering the matrix 
element correction does depend on more parameters^ than 
the one parameter of the shower kernel (corresponding to 
the dipole invariant mass). 



In a realistic application these are not two random num- 
bers needed for the Born process, but indeed six, since photon 
radiation is generated of each incoming lepton, requiring two 
random numbers per incoming lepton. 
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Parameter 


Range 


Description 


a.(Mi) 


0.1-0.13 


Input a s at Z mass. 


HIR.FF 


0.5 GeV - 2.0 GeV 


Infrared cutoff for final-final dipoles 


Msoft, FF 


0.0 GeV - 1.2 GeV 


Soft scale for final-final dipoles 


m g , c 


0.67 GeV - 3.0 GeV 


Gluon constituent mass 


Clmax 


0.5 GeV - 10 GeV 


Maximum cluster mass 


Clpow 


0.0-10.0 


Cluster mass exponent 


Clsmr 


0.0-10.0 


Cluster direction smearing 


Psplit 


0.0-1.4 


Cluster mass splitting parameter 



Table 1. The parameters varied for the fit to LEP data. 



Parameter 


LO 


NLO 


a s (Mi) 


0.113185 ±0.007281 


0.117550 ±0.005053 


fJ-IR,FF 


(1.416023 ±0.306430) GeV 


(1.245196 ±0.226821) GeV 


jUsoft, FF 


(0.242725 ± 0.202069) GeV 


0.0 GeV0 


m 9,c 


(1.080386 ±0.499546) GeV 


(1.007680 ±0.265565) GeV 


Clmax 


(4.170320 ± 0.589504) GeV 


(3.664004 ± 0.639504) GeV 


Clpow 


5.734681 ± 1.006965 


5.687022 ± 0.869322 


Clsmr 


4.548755 ± 2.350193 


3.115744 ± 2.436793 


-Psplit 


0.765173 ± 0.074008 


0.771329 ± 0.074248 


Table 2. Parameters for LO and NLO fits to LEP data. 


Parameter 


LO 


NLO 


I^IR,FI 


(0.796205 ± 0.333340) GeV 


(0.718418 ± 0.210448) GeV 


Msoft, FI 


(1.355894 ± 0.432515) GeV 


(1.003714 ±0.252398)GeV 


Table 3. Parameters for LO and NLO fits to HERA data. 


Parameter 


LO 


NLO 


I^IR.II 


0.367359 GeV 


0.275894 GeV 




0.205854 GeV 


0.254028 GeV 


-*4_1_ .valence 


1.68463 GeV 


1.26905 GeV 


.sea 


1.29001 GeV 


1.1613 GeV 



Table 4. Parameters for LO and NLO fits to the CDF Drell-Yan data. 



3 This parameter was predicted negative by Professor though consistent with zero and has thus been fixed. 



